function [] = b1t1qmtSimu_Stage1(outputDir,mt_measurement,vol_dims,vol_resolution,protocol, sled_trs, sled_excite_flips)
%B1T1QMTSIMU_STAGE1 Creates Minc files for the qMT fitting to be performed, according to the B1 and T1 maps defined within this file. 
%
%   -The qMT images were are flat (every pixel the same value) and equal
%   to the simulated data from main_qmt_simu.m.
%
%   -The B0 measurements are flat and equal to zero.
%
%   -The B1 map is a 1D gradient along the left-right direction.
%
%      B1_min           B1_max
%           _ _ _ _ _ _ _ _ 
%          |               |
%          |               |
%          |               |
%          |               |
%          |               |
%          |               |
%          |               |
%          |_ _ _ _ _ _ _ _|
%    
%
%   -The T1 map is a 1D gradient along the up-down direction.
%           _ _ _ _ _ _ _ _ 
%  T1_min  |               |
%          |               |
%          |               |
%          |               |
%          |               |
%          |               |
%          |               |
%  T1_max  |_ _ _ _ _ _ _ _|
%
%
%   -The dR1 images are required by the qMT fitting software, but are set 
%   to 0.
%
%   ***NOTE THAT NIAK REORIENTS THE IMAGES, SO THAT THE MINC FILES WILL NOT
%   HAVE THE GRADIENTS IN THE SAME DIRECTION***
%
%   Author: Mathieu Boudreau
%   Date Created: June 2014
%   Date Last Modified: June 20th 2014.
%



%% Create qMT Minc volumes
%
%Set No-MT volumes
qmt_m0_vol=ones(vol_dims);

if length(sled_trs)==1 %Use for UK protocol
    qMT_MINC_M0=minc;
    qMT_MINC_M0.fileName=[outputDir 'qMT/qmt_0.mnc'];
    qMT_MINC_M0=qMT_MINC_M0.setVolume(qmt_m0_vol);
    qMT_MINC_M0=qMT_MINC_M0.setResolution(vol_resolution);
    qMT_MINC_M0=qMT_MINC_M0.setTR(sled_trs);
    qMT_MINC_M0=qMT_MINC_M0.setFA(sled_excite_flips);
    qMT_MINC_M0.saveMinc;
else
    for ii=1:length(sled_trs)
    qMT_MINC_M0=minc;
    qMT_MINC_M0.fileName=[outputDir 'qMT/qmt_0_tr_' num2str(sled_trs(ii)) '.mnc'];
    qMT_MINC_M0=qMT_MINC_M0.setVolume(qmt_m0_vol);
    qMT_MINC_M0=qMT_MINC_M0.setResolution(vol_resolution);
    qMT_MINC_M0=qMT_MINC_M0.setTR(sled_trs(ii));
    qMT_MINC_M0=qMT_MINC_M0.setFA(sled_excite_flips(ii));
    qMT_MINC_M0.saveMinc;                
    end
end


% Set MT weighted volumes
for qmtCount=1:length(mt_measurement)
   qMT_MINC(qmtCount)=minc;
   qMT_MINC(qmtCount).fileName=[outputDir 'qMT/qmt_' num2str(qmtCount) '.mnc'];
   qMT_MINC(qmtCount)=qMT_MINC(1,qmtCount).setVolume(mt_measurement(qmtCount).*ones(vol_dims));
   qMT_MINC(qmtCount)=qMT_MINC(1,qmtCount).setResolution(vol_resolution);
   qMT_MINC(qmtCount)=qMT_MINC(1,qmtCount).setTR(protocol(qmtCount,4));
   qMT_MINC(qmtCount)=qMT_MINC(1,qmtCount).setFA(protocol(qmtCount,5));
   qMT_MINC(qmtCount).saveMinc;
end

%% Create B0 Minc volume
%

b0_vol=zeros(vol_dims);

B0 = minc;
B0.fileName=[outputDir 'B0/b0.mnc'];
B0=B0.setVolume(b0_vol);
B0=B0.setResolution(vol_resolution);
B0.saveMinc;

%% Create B1 Minc volume
%

b1_range=[0.5 2];
%b1_range=[1 1];

b1_vol=repmat(linspace(b1_range(1),b1_range(2),vol_dims(1)),vol_dims(2),1);

B1 = minc;
B1.fileName=[outputDir 'B1/b1.mnc'];
B1=B1.setVolume(b1_vol);
B1=B1.setResolution(vol_resolution);
B1.saveMinc;

%% Create T1 Minc volume
%

t1_range=[0.1 4];
%t1_range=[0.9 0.9];

t1_vol=repmat(linspace(t1_range(1),t1_range(2),vol_dims(1))',1,vol_dims(2));
T1 = minc;
T1.fileName=[outputDir 'T1/t1.mnc'];
T1=T1.setVolume(t1_vol);
T1=T1.setResolution(vol_resolution);
T1.saveMinc;

%% Create dR1 Minc volume
%

dr1_vol=zeros(vol_dims);

dR1 = minc;
dR1.fileName=[outputDir 'T1/dr1.mnc'];
dR1=dR1.setVolume(dr1_vol);
dR1=dR1.setResolution(vol_resolution);
dR1.saveMinc;

end

